#fit random intercepts and see how this affects \hat \rho

rm(list=ls())
library(epicalc)
zap()
#Spatial spec correct, omitted variable in X\betas (correlated with \rho)

sar.f <- function(es, x, y, w){
    rho <- es[1]
    s2 <- es[2]^2
    if (s2<0) return(-1e20)
    k <- ncol(x)
    b <- as.matrix(es[3:(k+2)])
    llik <- 0
    nd <- nrow(w)
    n <- nrow(x)
    d <- n/nd
    y <- as.matrix(y)
    
    for(i in 1:d){
        A <- diag(1, nd, nd)-rho*w[,,i]
        dA <- det(A)
        if (dA <=0){
            llik <-  -1e20 
            break
        } else {
            aux <- log(dA) -nd/2 * log(s2)- (1/(2*s2)) * t(A %*% y[((i-1)*nd+1):((i-1)*nd+nd),] - as.matrix(x[((i-1)*nd+1):((i-1)*nd+nd),])%*%b) %*% (A %*% y[((i-1)*nd+1):((i-1)*nd+nd),]- as.matrix(x[((i-1)*nd+1):((i-1)*nd+nd),])%*%b)
            llik <- llik + aux
        }
    }
    return(llik)
}

sar <- function(x, y, w, stval,vars=c(""), routit=100000){
    mod <- optim(stval, sar.f, hessian = TRUE, method = "BFGS", 
            control = list(fnscale = -1, REPORT = 10, 
            reltol = 1e-16, trace = 1, 
            maxit = routit), x=x, y=y, w=w)
    vacov <- as.matrix(solve(-1*mod$hessian, tol=1e-24))
    se <- as.matrix(sqrt(diag(vacov)))
    zstat <- as.matrix(mod$par / se)
    pvalues <- round(2*pnorm(abs(zstat),lower.tail=FALSE),4)
    
    cat("results in order: rho, s, b \n")
    cat("parameter estimates: ", mod$par, "\n")
    cat("standard error: ", se, "\n")
    cat("z-statistics: ", zstat, "\n")
    cat("p-values: ", pvalues, "\n")
    list(llik=mod$value, par=mod$par, se=se, zstat=zstat, 
        pvalues = pvalues, 
        hessian=mod$hessian, vacov=vacov
        )
}



genY <- function(nd, d, sig2, sigint, locint, rho, b0, b1, b2){ #Generate data from correct DGP

    wall.ev <- array(NA, c(nd, nd, d))
    wraw <- matrix(c(rep(1/(nd-1),nd)),nd, nd, byrow=TRUE)
    diag(wraw) <- 0
    wall.ev[,,1:d] <- wraw
        
    #draw random intercepts
    int <- rnorm(d,locint, sqrt(sigint))
        
    ytrue <- rep(NA, d*nd)
    inttrue <- rep(NA, d*nd)    #collect values of random intercept
    xother1 <- rep(NA, d*nd)        
    xother2 <- rep(NA, d*nd)        

    for (j in 1:d){
        iA <- solve(diag(1,nd)-rho * wall.ev[,,j])
        xotheraux1 <- rnorm(nd)
        xotheraux2 <- rnorm(nd)
        yaux <- iA %*% (int[d]+b0 + b1 * xotheraux1 + b2 * xotheraux2) + iA%*%rnorm(nd, sd=sqrt(sig2))
        ytrue[(1+(j-1)*nd):(nd+(j-1)*nd)]<- yaux
        inttrue[(1+(j-1)*nd):(nd+(j-1)*nd)] <- rep(int[j], nd)
        xother1[(1+(j-1)*nd):(nd+(j-1)*nd)]<- xotheraux1
        xother2[(1+(j-1)*nd):(nd+(j-1)*nd)]<- xotheraux2
    }
    return(list(ytrue=ytrue, inttrue=inttrue, xother1=xother1, xother2=xother2, wall.ev=wall.ev))
}


#repetitions
r <- 1000
    
#specify parameters
sig2 <- 1.0
sigint <- 2
locint <- 2
rho <- -.5
b0 <- .5
b1 <- 1
b2 <- -2

nd<-10
d<-100
set.seed(021175)

mod.true <- vector("list", r)
mod.lack <- vector("list", r)
i <- 1
while (i <= r){    
    mcdat <- genY(nd, d, sig2, sigint, locint, rho, b0, b1, b2)
    xfull <- cbind(1, mcdat$inttrue, mcdat$xother1, mcdat$xother2)
    xlack <- cbind(1, mcdat$xother1, mcdat$xother2)
    y <- mcdat$ytrue
    stvfull <- runif((ncol(xfull) + 2))*.1
    stvlack <- runif((ncol(xlack) + 2))*.1
    mod.true[[i]] <- try(sar(xfull, y, mcdat$wall.ev, stval=stvfull, vars=c("const", "ld", "other")), TRUE)
    ect <- 1 #initiate error counter
    while(class(mod.true[[i]]) == "try-error"){
            cat("\n\n restarting repetition ",i,", ",ect,". attempt \n\n")
            stvfull <- runif((ncol(xfull) + 2))*.1
            mod.true[[i]] <- try(sar(xfull, y, mcdat$wall.ev, stval=stvfull, vars=c("const", "ld", "other")), TRUE)
            ect <- ect+1
            if (ect == 10) {
                i <-  r+1
                cat("\n\n\n\n ROUTINE FAILED BC OF TRUE MODEL\n")
                break
            }
    }    
    mod.lack[[i]] <- try(sar(xlack, y, mcdat$wall.ev, stval=stvlack, vars=c("const", "other")), TRUE)
    ect <- 1 #initiate error counter
    while(class(mod.lack[[i]]) == "try-error"){
            cat("\n\n restarting repetition ",i,", ",ect,". attempt \n\n")
            stvlack <- runif((ncol(xlack) + 2))*.1
            mod.lack[[i]] <- try(sar(xlack, y, mcdat$wall.ev, stval=stvlack, vars=c("const", "other")), TRUE)
            ect <- ect+1
            if (ect == 10) {
                i <-  r+1
                cat("\n\n\n\n ROUTINE FAILED BC OF FALSE MODEL\n")
                break
            }
    }
    cat("\n\n finished repetition ", i, " \n\n")
    i <- i + 1
}

getpars <- function(i, model){
    model[[i]]$par
}
getllik <- function(i, model){
    model[[i]]$llik
}
pars.true <- unlist(sapply(1:1000, getpars, mod.true), recursive=FALSE)
pars.true <- matrix(pars.true, r, length(stvfull), byrow=TRUE)

pars.lack <- unlist(sapply(1:1000, getpars, mod.lack), recursive=FALSE)
pars.lack <- matrix(pars.lack, r, length(stvlack), byrow=TRUE)

llik.true <- unlist(sapply(1:r, getllik, mod.true), recursive=FALSE)
llik.lack <- unlist(sapply(1:r, getllik, mod.lack), recursive=FALSE)

quantile(pars.true[,1], c(.025,.5,.975))
quantile(pars.lack[,1], c(.025,.5,.975))

bias.true <- mean(pars.true[,1]-rho)
bias.lack <- mean(pars.lack[,1]-rho)

mean(llik.true)
mean(llik.lack)


#true model has better fit:
#ttest of likelihood
t.test(llik.true,llik.lack,paired=TRUE)

#llratio test of mean llik
dgf=1
r <- -2 * (mean(llik.lack)-mean(llik.true))
p <- 1-pchisq(r, dgf)
cat("degrees of freedom: ", dgf, "\n")
cat("likelihood ratio statistic: ", r,"\n")
cat("p value: ", p, " (under H0 that the restricted model is correct)\n\n")
# number of lliks greater for full model
check<- llik.true>llik.lack
sum(check)

setwd("C:/Users/Martin Steinwand/Documents/Project Lead Donor/Analysis/R&R.IO")
#save(mod.lack, file="mcld.omit.out")
#save(mod.true, file="mcld.full.out")
